#include<bits/stdc++.h>
using namespace std;
const int N = 1e7+10;
const int MOD = 1e9+7;
const long long MOD1 = 1ll*MOD*MOD;
const __int128 one = 1;
long long mu[N];
long long phi[N];
int minP[N];
vector<int> p;
int qpow(int x,int p){
	int res = 1;
	while(p){
		if(p & 1){
			res = 1ll*res*x%MOD;
		}
		x = 1ll*x*x%MOD;
		p >>= 1;
	}
	return res;
}
long long bigphi(long long vl){
	long long res = 1;
	for(int i = 0;1ll*p[i]*p[i] <= vl && i < p.size();i++){
		if(vl % p[i] == 0){
			res = res*(p[i] - 1);
			vl /= p[i];
			while(vl % p[i] == 0){
				vl /= p[i];
				res = res*p[i];
			}
		}
	}
	if(vl != 1){
		res = res*(vl - 1);
	}
	return res;
}
void init(int ma){
	mu[1] = 1;
	phi[1] = 1;
	for(int i = 2;i<=ma;i++){
		if(!minP[i]){
			minP[i] = i;
			p.push_back(i);
			mu[i] = -1;
			phi[i] = i - 1;
		}
		for(int j = 0;p[j] != minP[i] && 1ll*i*p[j] <= ma;j++){
			mu[i * p[j]] = -mu[i];
			minP[i * p[j]] = p[j];
			phi[i * p[j]] = phi[i]*(p[j] - 1);
		}
		if(1ll*minP[i]*i <= ma){
			minP[i * minP[i]] = minP[i];
			phi[i * minP[i]] = phi[i]*minP[i];
		}
	}
	for(int i = 1;i<=ma;i++){
		mu[i] = (mu[i] + MOD1 + mu[i - 1])%MOD1;
	}
}
typedef vector<vector<long long>> mat;
mat operator*(const mat &a,const mat &b){
	int n = a.size();
	mat res(n,vector<long long>(n,0ll));
	for(int i = 0;i < n;i++){
		for(int j = 0;j < n;j++){
			for(int k = 0;k < n;k++){
				res[i][j] = (res[i][j] + one*a[i][k]*b[k][j]%MOD1)%MOD1;
			}
		}
	}
	return res;
}
long long getf(const mat &m,long long n){
	n--;
	mat res(2,vector<long long>(2,0ll)),base = m;
	res[0][0] = res[1][1] = 1;
	while(n){
		if(n & 1){
			res = base*res;
		}
		base = base*base;
		n>>=1;
	}
	return res[0][1];
}
const long long i6 = 833333345000000041ll;
void solve(long long n,int a){
	long long r1 = 0;
	long long r2 = 0;
	int l = 1,r;
	for(;l <= a;l = r + 1){
		r = min(a/(a/l),a);
		r1 = (r1 + one*(mu[r] - mu[l - 1] + MOD1)%MOD1*(a/l)%MOD1*(a/l)%MOD1*(a/l)%MOD1)%MOD1;
		r2 = (r2 + one*(mu[r] - mu[l - 1] + MOD1)%MOD1*(a/l)%MOD1*(a/l)%MOD1)%MOD1;
	}
	r1 = one*(r1 + 2 + one*r2*3%MOD1)%MOD1*i6%MOD1;
	mat sf(2,vector<long long>(2));
	sf[0][0] = MOD1-1;
	sf[0][1] = one*r1*(r1 - 1)%MOD1;
	sf[1][0] = 0;
	sf[1][1] = (r1 - 1 + MOD1)%MOD1;
	long long ans = 0;
	long long pn = bigphi(n);
	for(int i = 1;1ll*i*i <= n;i++){
		if(n % i == 0){
			//cerr<<(pn/phi[i])*(phi[__gcd(n/i,(long long)i)]/__gcd(n/i,(long long)i))<<'\n';
			//cerr<<pn%phi[i]<<'\n';
			ans = (ans + one*pn*phi[__gcd(n/i,(long long)i)]/phi[i]/__gcd(n/i,(long long)i)*getf(sf,i)%MOD1)%MOD1;
			if(n / i != i){
				ans = (ans + one*phi[i]*getf(sf,n/i)%MOD1)%MOD1;
			}
		}
	}
	if(n % MOD){
		ans %=MOD;
		ans = ans*qpow(n%MOD,MOD - 2)%MOD;
	}else{
		ans /= MOD;
		ans %= MOD;
		ans = ans*qpow(n/MOD,MOD - 2)%MOD;
	}
	cout<<ans<<'\n';
}
int main(){
	int T;
	cin>>T;
	vector<pair<long long,int>> que;
	int ma = 0;
	while(T--){
		long long n;
		int a;
		cin>>n>>a;
		que.push_back(make_pair(n,a));
		ma = max(a,ma);
		ma = max(ma,(int)(sqrt(n)));
	}
	init(ma);
	for(pair<long long,int> qs : que){
		solve(qs.first,qs.second);
	}
}